#Example of reading .bin and .tif files, and making plots

#needed libraries (mostly for plotting)
library(ggplot2)
library(raster)
library(RColorBrewer)
library(grid)
library(gridExtra)

#Set working directory (update for your machine as needed)
setwd("~/Documents/ThesisPaper/FinalPaperSubmission/Hein_EtAl_Supplement")

#Load cell coordinates, set model dimensions
Coords = read.csv("./Cell_Coords.csv", header = 0)
colnames(Coords)=c("UTM_X","UTM_Y")
nx = 1200
ny = 1124

#Prepare data frame for input
Runs = c("Base1","Crops3","Hot3", "Dry3","HotDry3","Worst3")
nRuns = length(Runs)
ET = as.data.frame(matrix(nrow=1200*1124,ncol = nRuns))
colnames(ET)=Runs

#Read in the .tif file and reshape to a vector
Veg_matrix = as.matrix(raster("./Veg_Classes.tif"))
Veg_vector = rep(0,nx*ny)
p = 1
for(i in ny:1) #start at the bottom and work to the top
{
  for(j in 1:nx) #work right to left
  {
    Veg_vector[p]<-Veg_matrix[i,j]
    p=p+1
  }
}

#Load the .bin files
for(i in 1:nRuns)
{
  filename = paste("./ET/",Runs[i],"qflx_total.yearly.bin", sep="")
  to.read = file(filename, "rb")
  header=readBin(to.read, integer(), endian = "little", size=4, n=3)
  ET[,i]=readBin(to.read, double(), endian = "little", size=8, n=1200*1124) #The model has 1200 columns and 1124 rows
  close(to.read)
}

#make a plot of ET 
#note, we convert units to mm by multiplying by 1000
ETPlot = ggplot(data=Coords, aes(y=UTM_Y, x=UTM_X)) +
  geom_raster(aes(fill=ET$Base1*1000)) + 
  coord_equal() + theme_bw()+ labs(title="Baseline ET", y="", x="")+
  scale_fill_gradient("mm")+
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        axis.ticks = element_blank(), axis.text = element_blank() )
ETPlot

#make a plot of ET anomaly
#note, we convert units to mm by multiplying by 1000
ETAnomPlot = ggplot(data=Coords, aes(y=UTM_Y, x=UTM_X)) +
  geom_raster(aes(fill=(ET$Worst3-ET$Base1)*1000)) + 
  coord_equal() + theme_bw()+ labs(title="Worst Case ET Anomaly", y="", x="")+
  scale_fill_gradient2("mm")+
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        axis.ticks = element_blank(), axis.text = element_blank() )
ETAnomPlot

#make a plot of vegetation
Vegetation = as.factor(Veg_vector)
#    1      4      5      6      8      9     10     12     13     14     15     16     18 
#78170  21443   9145 214789   5908 123898 506398 227215   3841 152192     40    678   5083 
#Evergreen needleleaf	1--#deciduous broadleaf forests	4--#mixed forests	5
#closed shrublands	6--#woody savannas	8--#savannas	9--#grasslands	10
#croplands	12--#urban and built-up lands	13--#cropland / natural vegetation mosaics	14
#snow and ice	15--#barren or sparsely vegetated	16--#bare soil	18
my_labels3 = c("Evergreen needleleaf","Deciduous Broadleaf","Mixed Forests","Closed Shrublands","Woody Savannas","Savannas","Grasslands","Croplands","Urban","Crop Mosaics","Snow","Barren","Bare Soil")
cols = colorRampPalette(brewer.pal(8,'Dark2'))(13)
summary(Vegetation)
VegPlot = ggplot(data=Coords, aes(y=UTM_Y, x=UTM_X)) +
  geom_raster(aes(fill=Vegetation)) + 
  coord_equal() + theme_bw()+ labs(title="Vegetation Classes", y="", x="")+
  scale_fill_manual(values =cols, labels = my_labels3)+
  theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        axis.ticks = element_blank(), axis.text = element_blank() )
VegPlot
